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Abstract: In this paper, the Feynman path integral formulation of the continuous- 
continuous filtering problem, a fundamental problem of applied science, is investigated for 
the case when the noise in the signal and measurement model is additive. It is shown that 
it leads to an independent and self-contained analysis and solution of the problem. A con- 
sequence of this analysis is Feynman path integral formula for the conditional probability 
density that manifests the underlying physics of the problem. A corollary of the path in- 
tegral formula is the Yau algorithm that has been shown to be superior to all other known 
algorithms. The Feynman path integral formulation is shown to lead to practical and imple- 
mentable algorithms. In particular, the solution of the Yau PDE is reduced to one of function 
computation and integration. 
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1. Introduction 
1.1 Motivation 

The fundamental dynamical laws of physics, both classical and quantum mechanical, are de- 
scribed in terms of variables continuous in time. The continuous nature of the dynamical 
variables has been verified at all length scales probed so far, even though the relevant dy- 
namical variables, and the fundamental laws of physics, are very different in the microscopic 
and macroscopic realms. In practical situations, one often deals with macroscopic objects 
whose state variables are phenomenologically well-described by classical deterministic laws 
modified by external disturbances that can be modelled as random noise. It is therefore 
natural to consider the problem of the evolution of a state of a signal of interest described 
by a continuous-time stochastic dynamical model. Roughly speaking, the state of the sys- 
tem is described by a noisy version of a deterministic dynamical system termed the state 
model, i.e., the dynamics is governed by a system of first-order differential equations in the 
state variable (x(t)) with an additional contribution due to noise that is random(i>'(t)). The 
noise in the state model is referred to as the signal noise. If the noise is Gaussian (or more 
generally multiplicative Gaussian) the state process is a Markov process. Since the process 
is stochastic, the state process is completely characterized by a probability density function. 
The Fokker-Planck-Kolmogorov foward equation (FPKfe) describes the evolution of this prob- 
ability density function (or equivalently, the transition probability density function) and is 
the complete solution of the state evolution problem. For an excellent discussion of stochastic 
processes and filtering theory, see . 

However, in many applications the signal, or state variables, cannot be directly observed. 
Instead, what is measured is a nonlinearly related stochastic process (y(t)) called the mea- 
surement process. The measurement process can often be modelled as yet another continuous 
stochastic dynamical system called the measurement model. Loosely speaking, the observa- 
tions, or measurements, are discrete-time samples drawn from a different system of noisy first 
order differential equations. The noise in the measurement dynamical system is referred to 
as measurement noise. For simplicity, the signal and measurement noise are often assumed 
to be independent. 
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The problem of optimal nonlinear filtering is to estimate the state of a stochastic dy- 
namical system optimal under a certain criterion, given the observations that are samples 
of a related stochastic measurement process. The conditional probability density function 
of the state parameters given the observations is the complete solution of the filtering prob- 
lem since it contains all the probabilistic information about the state process that is in the 
measurements and in the initial condition. This is the Bayesian approach, i.e., the a priori 
initial data about the signal process contained in the initial probability distribution of the 
stsite{u{to,x)) is incorporated into the solution. Given the conditional probability density, 
optimality may be defined under various criterion. Usually, the conditional mean, which is 
the least mean-squares estimate, is studied due to its richness in results and mathematical 
elegance. The solution of the optimal nonlinear filtering problem is termed universal, if the 
initial distribution can be arbitrary. 



1.2 Fundamental Sochastic Filtering Results 

When the state and measurement processes are linear, the linear filtering problem was con- 
sidered by Kalman and Bucy|§, Q. The Kalman filter has been successfully applied to a 
large number of problems. 

Neverthless, the Kalman filter suffers from some major limitations. The Kalman filter is 
not optimal even for the linear filtering case if the initial distribution is not Gaussian^. In 
other words, the Kalman filter is not a universal optimal filter, even when the filtering problem 
is linear. In the general nonlinear case, the filter state is infinite dimensional, namely the 
whole conditional probability distribution function, and the Kalman filter cannot be optimal, 
although it may still work quite well if the nonlinearity is mild enough. 

The continuous-continous nonlinear filtering problem (i.e., continuous-time state and 
measurement stochastic processes) was studied in |, ||, @ and 0. This led to a stochastic 
differential equation, called the Kushner equation, for the conditional probability density 
in the continuous-continuous filtering problem. It was noted in|^], [^], and [|lO| that the 
stochastic differential equation satisfied by the unnormalized conditional probability density, 
called the Duncan-Mortensen-Zakai (DMZ) equation, is linear and hence considerably simpler 
than the Kushner equation. In and [^] were derived the robust DMZ equation, a partial 
differential equation (PDE) that follows from the DMZ equation via a gauge transformation. 

A disadvantage of the robust DMZ equation is that the coefficients depend on the mea- 
surements. Thus, one does not know the PDE to solve prior to the measurements. As a 
result, real-time solution is impossible. In [13|, it was proved that the robust DMZ equation 



is equivalent to a partial differential equation that is independent of the measurements, which 
is referred to as the Yau Equation (YYe) in this paper (see discussion at the end of Section 
for why it is being distinguished from the FPKfe). Specifically, the measurements only 
enter as initial condition at each measurement step. Thus, no on-line solution of a PDE 



^It may still be optimal for a linear system under certain criteria, such as minimum mean square error, but 
not a general criterion. 
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is needed; all PDE computations can be done off-line. As discussed later, there are other 
real-time solutions possible; however, the approximation error in those algorithms are not as 
well controlled. 

However, numerical solution of partial differential equations presents several challenges. 
A naive discretization may not be convergent, i.e., the approximation error may not vanish 
as the grid size is reduced. Alternatively, when the discretization spacing is decreased, it may 
tend to a different equation, i.e., be inconsistent. Furthermore, the numerical method may be 
unstable. Finally, since the solution of the YYe is a probability density, it must be positive 
which may not be guaranteed by the discretization. 



A different approach to solving the PDE was taken in |14] and [15]. An explicit expression 
for the fundamental solution of the YYe as an ordinary integral was derived. It was shown 
that the formal solution to the YYe may be written down as an ordinary, but somewhat 
complicated, multi-dimensional integral, with the integrand an infinite series. In addition, an 
estimate of the time needed for the solution to converge to the true solution was presented. 

1.3 Outline of the Paper 

In this paper, the Feynman path integral (FPI) formulation is employed to tackle the continuous- 
continuous nonlinear filtering problem. Specifically, phrasing the stochastic filtering problem 
in a language common in physics, the solution of the stochastic filtering problem is presented. 
In particular, no other result in filtering theory (such as DMZ equation, robust DMZ equation, 
etc) is used. The path integral formulation leads to a path integral formula for the transition 
probability density for the general aditive noise case. A corollary of the FPI formalism is 
the path integral formula for the fundamental solution of the YYe and the Yau algorithm — 
the central result of nonlinear filtering theory. It is noted that this paper provides detailed 
derivation of results that were used in In view of the wide applicability of the results 
derived here, this paper provides a brief summary of the central results of filtering theory and 
pedagogical discussion of Feynman path itegral methods that may be found to be useful by a 
wider audience, i.e., people specializing in different areas of applied science and engineering. 



The Feynman path integral was introduced by Feynman in quantum mechanics 1 17, 18] 



It has had a long and impressive history of results in theoretical physics. Inspired by an 



observation/remark of Dirac]19, 20], Feynman discovered the path integral representation of 
quantum physics ]17, 18]. The importance of path integrals grew (in the mid 1960's) when 
Fadeev and Popov used the path integrals to derive Feynman rules for the non-Abelian gauge 
theories. Its status was further enhanced when 't Hooft and Veltman used path integral 
methods to prove that the Yang-Mills theories were consistent (i.e., renormalizable). Since 
then, the path integral formulation of quantum field theory has led to extraordinary insights 
into quantum field theory, such as renormalization and renormalization group, anomalies, 
skyrmions, monopoles, instantons, and supersymmetry (see, for instance, the freely available 
text [^]). Much of the advance in the past 40 years in particle theory would not have been 
possible without the Feynman path integral. The path integral has also led to numerous 
insights into modern mathematics, especially in the work of Edward Witten. 
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A path integral formulation is especially attractive from a computational point of view. 
An excellent example is lattice quantum chromo dynamics (QCD) (see, for instance, p2| ). 
Specifically, lattice QCD has enabled high-precision calculations of observables in the non- 
Abelian gauge theory of strong interactions — an infinite dimensional, highly nonlinear prob- 
lem. The numerical results derived using the Feynman path integral based lattice QCD 
cannot be achieved by any other known technique. 

The following point needs to be emphasized to readers familiar with the discussion of 
standard filtering theory — Feynman path integral is different from the Feynman-Kac path 
integral. In filtering theory literature, it is the Feynman-Kac formalism that is used. The 
Feynman-Kac formulation is a rigorous formulation and has led to several rigorous results in 
filtering theory. However, in spite of considerable effort it has not been proven to be useful in 
the development of relaible practical algorithms with desirable numerical properties. It also 
obscures the physics of the problem. 

In contrast, it is shown that the Feynman path integral leads to formulas that are em- 
inently suitable for numerical implementation. It also provides a simple and clear physical 
picture. Finally, the theoretical insights provided by the Feynman path integral are highly 
valuable, as evidenced by numerous examples in modern theoretical physics (see, for instance, 
and shall be employed in subsequent papers. 

The outline of this paper is as follows. In Section |2|, aspects of continuous-continuous 
filtering are reviewed and some of the recent work on continuous-continuous filtering is also 
discussed. In the following section, some topics in path integrals are reviewed and the inter- 
pretation of transition probability density as ensemble averages. For more details on the path 
integral methods, see any modern text on quantum field theory, such as |21], and especially 
[p^ ] which discusses application of Feynman path integral to the study of stochastic pro- 
cesses. In Sections ^ and ^, the path integral formulas for the transition probability density 
are derived for the time-independent and time-dependent cases. In Section the partial dif- 
ferential equation satisfied by the path integral formula (derived for sampled measurements) 
in the previous sections is derived. It is shown that in the time-dependent case the partial 
differential equation is simply the YYe of continuous-continuous filtering. The conclusions 
are presented in Section 



2. Continuous-Continuous Filtering and the Yau Equation 



In this section, the main results of (continuous-continuous) nonlinear filtering theory are 
summarized. The most important results are the YYe and the Yau algorithm. In subsequent 
sections, it will be seen that they follow directly (and simply) from the Feynman path integral 
method. 
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2.1 The Continuous-Continuous Model 

The signal and observation model in continuous-continuous filtering is the following: 

(dx{t) = f{x{t),t)dt + e{x{t),t)dM{t), x{0)=xo, 

\dy{t) = h{x{t),t)dt + n{t)dw{t), y{0) = 0. 

Here x,v,y, and w are W\ M^, M™ and M'^ valued stochastic processes, respectively, and 
e{x(t),t) G MJ^^P and n{t) G R'"^'?. These are defined in the Ito sense. The Brownian 
processes v and w are assumed to be independent^ with p x p and q x q covariance matrices 
Q{t) and W{t), respectively. We denote n{t)W {t)ri^ (t) by R{t), a, m x m matrix. Also, / is 
referred to as the drift, e as the diffusion vielbein, and eQe^ as the diffusion matrix. 

In this paper, the additive noise case is considered, i.e., e{x{t),t) = e{t). Then, the Ito 
and Stratanovich forms are the same. In a future paper, the more general, multiplicative and 
correlated noise case shall be tackled. 

In this section, some of the relevant work on continuous-continuous filtering is summa- 
rized. Hence, it is assumed that n = p,m = q, f and h are vector-valued C°° smooth 
functions, e(x, t) is an orthogonal matrix-valued C°° smooth function, Q{t) is a n x n iden- 
tity matrix, and n{t) and R{t) are m x m identity matrices. No explicit time dependence is 
assumed in the model. 

2.2 The DMZ Stochastic Differential Equation 

The unnormalized conditional probability density, a(t,x) of the state given the observations 
{Y{s) : < s <t} satisfies the DMZ equation: 

m 

da{t,x) = ^Qa{t,x)dt + ^^ia{t,x)dyi{t), a{^,x) = uq. (2.2) 



Here 



i=l 



i=l ' i=l 1=1 1=1 

d 1 A 9' 



--E£-(/<wo4Eik-5E'^(^ 

i=l i=l ' i=l 



where is the zero-degree differential operator of multiplication by hi{x), i = 1, . . . ,m, ao 
is the probability density of the initial time to, and 

A ^ ^ - fi{x), (2.4) 



1=1 i=l i=l 



^Hence, the large body of work on the correlated noise case is not considered in this paper 
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The DMZ equation is to be interpreted in the Stratanovich sense. Note that 

Di = ^-U{x), (2.5) 

OXi 



where = F,(x) = / fi{t)dt. 
OXi Jo 



Hence, 



,2 



A^ = e^»|,e-^% (2.6) 



and 



^o = ^Ee^»^e---.(-)j. (2.7) 
2.3 The Robust DMZ Partial Differential Equation 



Following [^] and |12] introduce a new unnormalized density 

u{t, x) = exp ^- hi{x)yi{t)^ a{t, x). (2.8) 
Under this transformation, the DMZ SDE is transformed into the following time- varying PDE 

i=l ^ i=l \ j=l I 



( X 1 '''' 1 '''' ''' c\i^ 

E l^(-) + ^E^^(-) - ^E^^(*)^^^(-) +EE^^w/.(-)£'(^ 
i=l * 1=1 i=l 1=1 j=l ^ 



2 

i,j=i k=i 

^ u{0,x) =aoix). 

(2.9) 

This is called the robust DMZ equation. Here A is the Laplacian. The solution of a PDE 
when the initial condition is a delta function is called the fundamental solution. 

2.4 Universal Linear and Finite Dimensional Filtering I: Lie Algebra Method and 
the Maximal Rank Case 

The Kalman filter is not a universal optimal linear filter since it assumes that the initial 
distribution is Gaussian, although it may still be optimal under certain criteria such as minu- 
mum variance. In the general case, the DMZ equation needs to be considered. The use of Lie 
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algebraic methods to solve the DMZ equation had been proposed in [^] and |2^. The Lie 
algebraic viewpoint has been useful in giving a deeper understanding of the DMZ equation. 

The estimation algebra E associated with the continuous-continuous model is the Lie 
algebra generated by ^o,^i{x), . . . ,^m{x). It is said to be of maximal rank if for any 
1 < i < n there exist constants Cj such that Xi + Ci are in E. When m > n, if h{x) = Hx, 
where H is of maximal rank, then E is of maximal rank. Note that the explicit recursive 
solution for the Benes filtering system (see |26|]) was originally derived only for the maximal 
rank case. In some cases, maximal rank condition is assumed in deriving the explicit solutions. 

If the linear filtering system is completely reachable, controllable, and observable, the 
estimation algebra is 2n + 2 dimensional with basis 

d d 

J^'O,— ,Xu...,Xn,l. (2.10) 

The Wei-Norman approach to solve the robust DMZ equation, a time- variant PDE, is to note 
that the solution is of the form| 



n(t, x) = e^We''"W^" • • • g'-iW^ie^-W^" • • • e^iW^ie*^On(0, x), (2.11) 

where T{t),ri{t), . . . , r„(t), si(t), . . . ,Sn{t) satisfy a system of 2n -|- 1 ODEs. Therefore, the 
Wei-Norman approach requires the solution of a system of (2n -|- 1) ODEs, n linear PDEs, 
and an FPKfe type of equation. However, if the linear system is not completely reachable 
or completely observable, the basis of the estimation algebra is not known beforehand and 
therefore needs to to calculated. Hence, one cannot even write down the finite system of 
ODEs explicitly. 

2.5 Universal Linear and Finite Dimensional Filtering II: Direct Solution 

A direct method for solving the universal linear filtering problem was first presented in p8[| . 
It was shown that by a combination of a gauge transformation and a time-dependent spatial 
translation of the unnormalized conditional probability density, the robust DMZ equation for 
the linear filtering problem reduces to an on-line system of ODEs and an off-line FPKfe. This 
solution does not require controllability or observability assumptions, and is also universal as 
the initial distribution could be arbitrary. This solution was further simplified and extended 



in |29|1 and |30|]. 



The direct solution for the Benes case was also derived by noting that if / is a vector 
field of a potential function cj), then 



1 / " «2 \ 

^ E-'7w2-"'-^(-) • (2-12) 



If S,{t,x) = e '^a{t,x), then the DMZ equation is 



^(t,x) = ^(A-^(x))e(t,x) + X;^.e(i,x)^(t), (2.13) 

i=l 
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where A is the Laplacian operator. If 77 is a quadratic polynomial in x, the semigroup 
generated by the differential operator A — rj{x) is known and can be used to explicitly derive 
solutions to the equation when hi are linear in x. This technique is also related to the 



concept of equivalence of parabolic equations (see |28] for details). The solution was further 
generalized and simplified to the case of a more general finite-dimensional filter (Yau filter^ 
with quadratic r][x)). In summary, they showed that if one considers u{t,x) defined as 



u{t, x) 



exp 



c{t) + ai{t)xi - F{x + b{t)) 



u{t,x + b{t)), 



(2.14) 



then the robust DMZ equation for u{t) is reduced to a system of ODEs for a{t), b{t), c{t), and a 
FPKfe-type of equation for u{t, x). It has since been shown that even the off-line FPKfe-type 
of equation can be reduced to a system of linear ODEs [R^, 32]. Note that, like the Lie algebra 



solution, the state and measurement models are assumed to be explicitly time-independent. 
2.6 The Yau Equation 

The finite dimensional filtering problem is a very restrictive class of the general nonlinear 
filtering problem. Recently, it was proved that the real-time solution of the general nonlinear 
filtering problem can be obtained reliably ||l^, [^]. Let = {tq < n < • • • < = t} be 
a partition of the time interval [To,r], and let the norm of the partition be defined as 
\^k\ = sup;^<j<;, {ki ~ ''"i-ill- If '^li'tjx) satisfies the equation 



dui 
dt 



i=i « 

n 



1=1 

^ m 



dxi 



dui 
dxi 



{t,x) 



dhi 



1=1 



dx. 



i=l 



i=l 



i=l j=l 



dxj 



I E Ey*(^')2/j(ri)^(x)|^(x) ]ui{t,x), 



i,j=l k=l 
^ ui{ti_i,x) =Ui_i{ti^i,x), 



in the time interval r/_i ^ti, then the function ui{t,x) defined as 

ui{t,x) = exp (^Y^yi{Ti)hi{x)j ui{t,x) 



(2.15) 



(2.16) 



Yau filter is one where the drift in the signal model is the sum of a linear term and a gradient of 
a function 4'{x) and the measurement model is linear. The function 4>{x) is assumed to have no quadratic 
form piece; else it can be absorbed into the linear term. When the gradient term vanishes, the Kalman filter 
results. When the linear term vanishes and the gradient term is such that i] is quadratic, the Benes filter |26[ 
is obtained. 
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satisfies the parabolic partial differential equation 



i=l * i=l \i=l 1=1 / 

1=1 * 



(2.17) 



in tlie same time interval. The converse of the statement is also true. In [p4| , it was also 
noted that it is sufficient to use the previous observation, i.e., ui{t,x) satisfies Equation |2.15 
if and only if ui{t,x) defined as 

ui{t, x) = exp yi{Ti^i)hi{x)^ ui{t, x) (2.18) 

satisfies Equation 2.17 in the time interval r^-i < t < ti. We refer to Equations 2.16| ( p. 18 ) 



and Equation |2.17] as the post-measurement (pre- measurement) forms of the YYe. 



Observe that Equation 2.15 is obtained by setting y{t) to y(T;) in Equation 2^*. It was 
proved that the solution of Equation 2.15| approximates very well the solution of the robust 
DMZ equation (Equation ^l9|), i.e., it converges to u{t,x) in both pointwise sense and 



sense. Thus, solving Equation 2.£ is equivalent to solving Equation 2.17. Finally. 



C7(r,x) = ^^im ^Uk{n,x). 



(2.19) 



Thus, the solution of the YYe (as \^k\ ~^ 0) is the desired unnormalized conditional proba- 
bility density. 

Observe that when h{x) = 0, it is simply the FPKfe. However, unlike the FPKfe, the YYe 
does not satisfy the current conservation condition, i.e., the right-hand term is not a total 
divergence. This means that it does not conserve probability. This fundamental difference is 
traced to the fact that the FPKfe evolves the normalized probability density (and preserves 
the normalization), while the YYe evolves the unnormalized conditional probability density. 
Therefore, this distinction is made between the two equations in this paper. 

2.7 The Yau Algorithm 

We may summarize the real-time algorithm, based on both the pre- and post-measurement 
forms of the YYe, of Yau as follows. Suppose measurements are available at times 



<To <Tl <T2 < 



< Tk 



T. 



(2.20) 



■^In contrast, note that u{t,x) defined for solving tfie finite-dimensional Yau filter in Equation ^.14| in the 
aforementioned papers is equivalent to solving the robust DMZ equation for u{t, x) at each time instant. 
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We seek the solution Ui(t, x), which is the solution of the robust DMZ equation. Let the initial 
distribution be u{tq,x) = (To(x). Then, evolve the initial distribution to the first measurement 
instant, ti, using the YYe: 



ui{to,x) 



exp ( Z]'ili[yi('^o) - yj{T^i)]hj{x) ) cro(a;) (Pre-Measurement) 



6xp ( ~ yi('^o)]^j(2;) ) o-o(a;) (Post-Measurement). 

(2.21) 



The solution of equation 2.21 at time ri is ui{ti,x). Note that ui{ti,x) is given by 



, . \ewi-YALiyiiTo)hiix))ui{Ti,x) (Pre-Measuremnt) 
ni(ri,x) = < ^ (2.22) 

[exp (— ^™ X yi{'Ti)hi{x)) ui{ti,x) (Post-Measurement). 



Next, solve the YYe to the next measurement instant T2 with initial condition U2iTi,x), i.e., 

i=l * i=l * \i=l * i=l / 



■U2(ri,x) 



exp f EJli (yj(n) - yj{To)) hj{x)] ui{ti,x) (Pre-Measurement) 
( X^jli (yj (''"2 ) ~ Vj (''"1 ) ) (^) ) '^i{ti,x) (Post-Measurement) . 



(2.23) 



to obtain U2{t2,x). In fact, for i > 2, Ui{Ti,x) can be computed from Ui{Ti,x), where Ui{t,x) 
satisfies the equation 



1=1 * j=l Yj=1 ^ j=l 

exp (EJIi (yi('^i-i) - yjiTi-2)) hjix)^ ni_i(ri_i, x) (Pre-Measurement) 
exp (l]j=i (yj('^i) - yi(T"j-i)) hj{x)j ni_i(ri_i,x) (Post-Measurement). 



■Uj(ri_i,x) 



(2.24) 
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The initial condition in Equation 2.24 follows from noting that (since Mj(rj_i, x) = Ui-i{Ti^i, x)) 

/m \ 

Ui{Ti-i,x) = Ui{Ti-i,x)ex.p ^^yj{Ti)hj{x) , (2.25) 

(m \ / m \ 

— yj(rj_i)/ij (x) I ttj_i(Tj_i, x) exp yj(rj)/ij(x) I , (Pre-Measurement) 
j=i J \j=i J 

Ui{Ti-i,x) = nj(ri_i,a;)exp yj(Ti)/ij(2;)^ , (2.26) 

= exp — ^ yj(Tj_i)/ij(x) j nj_i(rj_i, x) exp '^^yj{Ti)hj{x) I (Post-Measurement). 
Thus the solution of the robust DMZ equation at the ith step is given by 



( ~ Yl'TLi yj{Ti-i)hji^) I Ui{Ti,x) (Pre-Measurement) 



u^in,x) = { } _l . _^ : . . , (2-27) 



6xp ( — YlT=i yj{Ti)^i{^) ) ^''^(''"ji (Post-Measurement). 



Note that the Yau algorithm is a recursive algorithm as it does not need any previous 
observation data. Furthermore, the YYe is independent of data and so can be computed 
off-line, and that the YYe is much simpler than the robust DMZ equation. Finally, note that 
the output of the Yau algorithm is the desired unnormalized conditional probability density. 

2.8 Other Real-Time Algorithms 

It is noted that there have been several other numerical methods proposed for solving the 



DMZ equation such as |3g, |37|, Some numerical examples were presented in |39(| , 



As noted in [^], most papers on the approximate solution of the DMZ equation deal with 
algorithms that are not (or, at least not easily) implementable. The most well-known of thsese 
"implementable" methods for the general filtering problem (rather than those applicable to 
only problems like finite-dimensional filters) are now briefly reviewed. 

Recall that the DMZ SDE for the unnormalized conditional probability density is given 

by 

n 

da{t, x) = Sfia{t, x)+Y^ ^i(T{t, x)dyi{t), (2.28) 

i=l 

where ^/ = + ^ hf{x), i.e., is the adjoint of the infinitesimal generator of the 
state Markov process, or the forward Fokker-Planck-Kolmogorov (FPK) operator. 
The simplest scheme is the explicit one-step Euler scheme: 

n 

a{t + A,x) = {1 + AJ(',)a{t,x) + ^hi{x)iyi{t + A) - yi{t))a{t,x), (2.29) 

i=l 
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Some other algorithms follow from observing that Equation 2.28 can be considered as the 
combination of a linear stochastic equation and the FPKfe: 



dai{t,x) = ai{t,x)'^hi{x)dyi(t), 

k=l 

da2{t, x) = ^/(J2(t, x)dt. 



(2.30) 



The standard idea of operator splitting for solving deterministic higher dimensional PDEs 
can then be applied to this case. Since the linear stochastic equation can be solved explicitly, 
the following two splitting-up approximations result: 



a{t + A,x) = < 



' / m 1 \ 

Tacxp I + A,x) - yi{t,x)\hi{x) - A- ^/i^(x) j (T{t,x), 

\i=l i=l J 

(m 1 ™ \ 
J2[y{t + A) - y{t)]h{x) hfix) TAa{t, x). 

i=l i=l / 



(2.31) 



Here Ta is the operator corresponding to the evolution via the FPK forward equation (FPKfe). 

The spectral separation scheme utilizes the fact that the conditional probability density 
can be written as a two-tensor basis with one of them obeying an observation-independent 
FPKfe type of equation (deterministic Hermite-Fourier coefficients) and the other a Wick 
polynomial that are completely defined by the observation process (see fs^ for details). This 
separation of parameters and observations also enables real-time implementation. It can be 
shown that it includes the explicit Euler one-step and splitting-up approximations. 

Finally, it will be seen that the path integral formulas derived here can be used to solve 
the FPKfe arising in these techniques. In fact, this formula already arises in the investigation 
of the continuous-discrete filtering problem (see |^0|). Note that in the Yau algorithm, the 
"prediction" part is accomplished via the YYe. In contrast, in the Splitting-up method, it is 
done using the FPKfe, as in continuous-discrete filtering theory. 

2.9 Why the Yau Algorithm is the most suitable algorithm 

The following two reasons explain why the Yau algorithm is superior to all other known 
algorithms: 

Firstly, recall that the coefficients of the robust DMZ equation (Equation 2.9) contain 
the measurements. This means that the partial differential equation to solve for continuous- 
continuous filtering is not even known prior to the measurements. In other words, the robust 
DMZ equation has to be solved on-line; its solution cannot be pre-computed. In contrast, in 
the YYe (Equation |2. 17 ), measurements are absent from the partial differential equation. The 
measurements only enter the initial condition at each measurement step. This implies that the 
YYe can be be solved off-line; there is no need for on-line solution of PDEs^. This makes real- 
time solution feasible even if state dimension is large. Specifically, suppose the measurements 



^This assumes that the measurement steps are equidistant or known. 
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are performed at equidistant time intervals. The set of all possible initial conditions may be 
expressed as a linear combination of basis functions. Then, one may precompute and store 
the solution to the YYe for the basis functions. Upon arrival of measurements, the initial 
condition incorporating the measurements is computed. This is projected onto the basis 
functions to ascertain the components. The solution (due to linearity of the PDE), till the 
next measurement, is then simply the sum of the computed off-line solutions for the basis 
functions. 

Secondly, all other known algorithms assume the signal and measurement model drift were 
bounded. In fact, error bounds derived for these methods depended (often exponentially) on 
the bounds so that there are severe time-step restrictions. Note that these results cannot 
even cover the Kalman filter (i.e., linear drift) case, which is a major limitation in solving 
practical problems. In contrast, the results of |Q and assume only that the drifts are 
nonlinear with linear growth. The errors are independent of the bounds on the signal or 
measurement model drift. In [^], this is further weakened: is is sufficient that the growth of 
the observation drift be greater than the growth of the signal model drift. 

It will be seen that the Feynman path integral approach naturally leads to the Yau 
algorithm. In that sense, the Feynman path integral formulation gives a physical explanation 
why the Yau algorithm is the natural choice. In addition, the path integral formula explicitly 
solves the PDE. 

2.10 Integral Formula for the Solution of the YYe 

A solution of the YYe in terms of ordinary integrals was presented in , jl^ . The formal 



solution derived is an infinite series in t (for details see the reference |15|). In addition, 
an estimate of the time interval on which this solution converges was presented. For the 
purposes of this paper, it is sufficient to note that it has better numerical properties (than 
solving PDEs) since integration can be computed reliably, as opposed to partial differential 
equations. However, the integrand is rather complicated (e.g., convolutions) and it does not 
appear to be easily implementable for a general model. We shall see that the path integral 
formula is considerably simpler to evaluate. 

3. Preliminaries 

In this section, a brief review of relevant aspects of functional methods common in theoretical 
physics is presented. For more details on these topics, the reader is referred to [p^ . 

3.1 A Brief Review of Functional Calculus 

The infinite dimensional generalization of a function is called a functional, F[x{t)], which 
gives a number for each function x{t). Loosely speaking, the functional derivative may be 
defined analogous to the ordinary derivative as follows: 

^^^""^^^^ = lim ^^^^^^ + '^^^ ~ ^'^^ ~ ^^"^^^^^ (31) 
Sx{t') e^o e ■ ^ ' ' 
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Alternatively, the functional derivative, s^^i^it)], of the functional with respect to 



variation of the function x{t) at s is defined by the following equation 
F[x{t) + m = F[x{t)] + I ^^[x{sMs')ds' + ■ 



(3.2) 



Note that the functional derivative of a functional is also a functional. Functional differenti- 
ation obeys the standard algebraic rules (linearity and Leibnitz's rule): 



Sx{t') 



(3.3) 



^6x{t')' 
.1=1 J 1=1 

^ [FMt)]F2[xm = FMt)]-^f2[x{t)] + F2[x{t)]^^FMt^ (3.4) 

Just as the derivative of x with respect to a; is 1, the functional derivative of x{t) with respect 
to x(t') is the unit matrix in infinite dimensions, namely the delta function: 

dx{t) 



5x{t') 



6{t - 1'). 



(3.5) 



The functional derivatives of a formal power series in function x{t) can be seen to be analogous 
to the ordinary derivative of a power series in the variable t. Specifically, it is straightforward 
to verify that if 



oo ^ 

F[<t)]=Y.-\ / dti---dtnF^''\tut2.---,tn)x{ti] 
n=0 



■x{tr. 



then 



FW(ti,t2,...,tr 



F 



(3.6) 



(3.7) 



X = Q 



In particular, a functional representable as a power series in functions, like the exponential 
functional, can be functionally differentiated using this result. Thus 

5 



6J{x] 



exp l^y" dx'J{x')A{x')^ = A{x) exp (^J dx'J{x')A{x')^ . (3.8) 



The functional delta function may be written as 



<5(/(x(t)) = J [dX{t)]exp (^i I dtXit)fixit)) 



(3.9) 



In the application at hand, determinants of operators (matrices) arise in Gaussian integra- 
tion in the infinite (finite) dimesional case. In particular, the operator det{d{x — y) + K{x, y)) 
needs to be evaluated for the problem at hand, for some operator K{x,y). It is assumed that 
traces of all powers of K exist. Prom the identity 



In det M = trlnM, 



(3.10) 
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it follows that 



lndet[l + i^]= J dxK{x,x) — - j dxidx2K{xi,X2)K{x2,x^) + ■ ■ ■ + (3-11' 



(_l)n+l 



n 



dXi ■ ■ ■ dXnK{xi,X2)K{x2,X3) ■ ■ ■ K{Xn, Xi) + 



Note that if K{x, y) = 6{x — y)K{x, y), only the first term is nonzero. 

Finally, recall the the sequence of identities from multivariable calculus 

{d"/(x)}5"(/(x)), 
{d^x}J6^{fix)), J = det(^^{x) 



Then, it follows that 



— Xc) 



J{x) ' 
5^{x-Xc)=5^{f{x))J{x). 



and 



a{x) = y {d"x}5"(/(x))J(x)a(x). 
Using the Fourier integral representation of the delta function, it is clear that 

{(i"x}{d"A}e^^-^(^)j(x)a(x). 



1 



a(x] 



(2vr)^ 



The functional analogue of Equation |3.12 

1 = 



IS 



^/(2;(t'))) det 
Finally, the Gaussian integral identity 



5f 



6x{t 



6{f{xm. 



(3.12) 



(3.13) 



(3.14) 



(3.15) 



(3.16) 



/ [d^x] exp I ^ Xi{A ^)ijXj - ^bjXj j = exp | ^ ^ biAijbj j , [d"x] = 



Y^(27r)"det^ 
(3.17) 



becomes 



-\tl J ^*^*(*) + dMt)Xi{t)j =exp U |; 1 6,(t)^,,(t)6,(t)j . 



(3.18) 
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3.2 Langevin Equation and Gaussian Measure 

In the path integral formulation, as in physics literature, the continuous-continuous model is 
interpreted as 



rx(t) =/(x(t),0 + e(x(t),t)i/(t), 
\y(i) =h{x{t),t)+n{t)n{t). 

Such equations are called Langevin equations. Heuristically, 

Gaussian noise process is represented by a path integral measure 



[dpiiym = [&v{t)] exp 



a,b=l 



5X1 



where is a real vector for each t. 

It is straightforward to verify that the mean is zero 

= 0, 

and the variance is Q{t)^: 



a,6=l 



(3.19) 



(3.20) 



(3.21) 



(3.22) 



(3.24) 



^a{t){Q-\t))am{t) + haUa{t) 



a, 6=1 



a=l 



Sbc{t) 6hd{t) 
= Qcd{mt-t'), 



exp I J dt' Y ba{t')Qab{t'Mt') 



a,b=l 



ba=0 



®This is also consistent with the previous discussion in terms of Ito calculus, i.e., dva{t)di'bit') = Qab{t)5tt'dt 
implies that 



(3.23) 



ba=0 
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and higher-order moments can be easily written down using Wick's theorem for bosonic 
(commuting) variables. The results are in accord with the expectation that the measure 
represents a Gaussian process with these first two moments. Finally, implicit in the use of 
these formal technques is the use of the Feynman convention, or symmetric discretization for 
the drift. 

3.3 Transition Probability Density and Signal and Measurement Ensemble Aver- 
ages 

Consider an ensemble of systems with state variables described by the Langevin equation. 
Due to the random noise, each system leads to a different vector x{t) that depends on time. 
Although only one realization of the stochastic process is ever observed, it is meaningful to 
speak about an ensemble average. For fixed times t = ti,i = 1,2, ... ,r, the probability density 
of finding the random vector x{t) in the (n— dimensional) interval Xi < x{ti) < Xi + dxi{l < 
i < r) is given by 

Wr{tr,xr;--- ;ii,xi) = ^[] '^"(x(ii) - x,) ^ , (3.25) 

where Xi is an n— dimensional column vector and (•) denotes averaging with respect to the 
signal model noise i'{t). The complete information on the random vector x(t) is contained 
in the infinite hierarchy of such probability densities. The quantity of interest here is the 
conditional probability density 

p{tr,Xr\tr-l,Xr-l, . . .]ti,Xi) = ((5"(x(i^) - x{tr))) U(t,_i)=x,_i,...,x(ti)=a;i > x{tr) = Xr 

(3.26) 

_ Wr{tr,Xr; . . ■ ;ti,Xi) 

/ WniU, Xr-,...; tl,Xl)d"-Xr ' 

Now the process described by the Langevin equation with 5— correlated Langevin force 
is a Markov process, i.e., the conditional probability density depends only on the value at the 
immediate previous time: 

p(tri Xr\tr—l, Xf—i, . . . ,ti, Xi) — pitni Xn\tn—li Xn—l) ■ (3.27) 

Hence the complete information for a Markov process is contained in the transition probabil- 
ities 

|, X W2{t2,X2;ti,Xi) 

p{t2, X2\h,xi) = , , (3.28) 

yyi[ti,xi) 

= (<5"(x(t2) - X2)) |x(ii)=xi- 

which satisfies the Chapman-Kolmogorov semigroup equation: 

P{h,Xz\ti,Xi) = j p{h,X2,\t2,X2)p{t2,X2\tl,Xi)dX2. (3.29) 
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It can be shown that the transition probabihty satisfies the Fokker-Planck-Kolmogorov 



forward equation (FPKfe) (see for instance, [42|) 



d "9 1 ^ 9^ 

^(t, x) = - — [fi{x{t),t)p{t, x)] + -Y,Y. Mxit),t)Qit))ijejkixit),t)p{t, x))] , 

i=l * i,k=l j=l * 

(3.30) 

= ^p{t,x). 

In pO|], the path integral formula for the fundamental solution for the FPkfe is derived and 
applied to the continuous-discrete filtering problem with additive (state model) noise. 

In continuous-continuous filtering, the measurement stochastic process needs to be in- 
corporated as well. Consider another ensemble of systems with state variables whose time 
evolution is governed by the measurement process. The measurement noise means that each 
system in the ensemble leads to a different time-dependent vector y{t). Thus, even though 
only one realization of the measurement stochastic process is observed, it is still meaningful 
to talk about an ensemble average of the measurement process (in addition to one over the 
state process). Thus, the quantity of interest in continuous-continuous filtering is 

P{t2,X2;y2\tl,Xi;yi) = (((5"(x(t2) -2;2)<J"'(y(t2) -y2)))^U(ti)=xi,j/(ti)=yi, (3-31) 

where (•)^ denotes averaging with respect to the measurement noise n{t). In the following 
sections, the path integral formulas for P{t2,X2;y2\ti,xi;yi) are derived. It shall be shown 
that the YYe plays the same role here that the FPKfe does in continuous-discrete filtering. 

Note that the transition probability density is the complete solution to the continuous- 
continuous filtering problem, since if the initial distribution is u(tj_i, x'|yi_i), then the evolved 
conditional probability distribution is 

u{t,x\Yi) = j P{t,x]y,\U_i,x';yi^i)u{U^i,x'\Y,_i) {(Tx'} . (3.32) 

In Section ^ it shall be shown that it leads to the Yau algorithm. 

4. Path Integral Formula: Implicit Time Dependence 
4.1 General Result 

In this section, the following signal and measurement models are considered: 




f{x{t))dt + du{t), X{Q) = Xo, Q{t) = Mnxn, 

h{x{t))dt + dti{t), y{0) = 0, R{t) = h^Imxm- 



(4.1) 



Here h,^ and are positive real numbers and Inxn is the nxn identity matrix. This is the 
same as the model in except in those references h^^ = hy = \, and an orthogonal vielbein 
matrix multiplied the state Brownian process. 
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Consider the member of the state ensemble with initial state x{to) = xq and final state 
x{t) = X and with signal noise sample history i'i{t),i = l,...,n. Then, from the state 
equation 



Xi{t) = fi{x{t)) + Vi{t), x{to) = Xq, x{t) = X. 



(4.2) 



Similarly, consider a member of the measurement ensemble with measurement noise sample 
history /Xj(f), i = 1, . . . , m. Then, the measurements satisfy the equation 



yk{t) = hk{x{t)) + fj,k{t), y{to) = yo, y{t) = y, 



(4.3) 



and k = 1, . . . ,m. The transition probability density is computed by averaging over the signal 
and measurement ensembles and imposing the condition that the signal and measurement 
Langevin equations are satisfied. That is, 

P{t,x;y\to,xo;yQ) = J [dp{iy{t))][dp{fi{t))] (4.4) 

X [dx(t)]JS''[x{t) - f{x{t)) - u{t)] 

x[dmJyS'^[m-h{x{t))-fi{t)] 

X ,5"[x(t) - x]^to)=.o^'^[yit) - y]\y{to)=yo- 

Here J and Jy are Jacobians that are computed below. Note that the measurement Langevin 
equation is imposed in the path integral expression. 

Based on the assumptions of the signal and measurement noise processes, it is clear that 



[dp{um = [^i.(t)]exp (^-±-j2iyiit)dt^ , 



[dp{ii{t))] = [&ii{t)] exp 



1 " /"* 



i4{t)dt 



(4.5) 
(4.6) 



Next, the Jacobians, J and Jy are computed. Observe that 



5xi{t') 



-[xi{t)-h{x{t))-v^{t)] 



d_ 



Sit-t'), 



8,^5{t-t')-e{t-t')^{x{t)) 



(4.7) 



smce 



±.e(t-t') = -5{t-t'), and ^^ = 5{t-t'). 



5x{t') 



(4.8) 
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The Jacobian J is thus given by 



5,^6{t-t')-9{t-t')^{x{t)) 



(4.9) 



det (^-^) det (^6ijd{t - t') - 9{t - t')^{x{t)) 
^det l6ij5{t-t')-e{t-t')^{x{t)) \ , 



where ^ is an irrelevant constant (independent of x) that may be absorbed in the measure. 

(4.10) 



From the identity det^ = exp(trln^), and using the identity in Equation 3.11 

In J = In det 



6,^6{t-t')-e{t-t')^{x{t) 



1=1 



The choice 6{0) = ^ ensures that commutativity of averaging and time differentiation ||23|| . 
Similarly, it is straightforward to see that the Jacobian Jy is trivial since y{t) is y— independent, 
i.e., 



- [Ut) - K{x{t)) - lli{t)] 



5yj{t 

and can be absorbed into the measure. 

Combining these results, it follows that 

ry(t)=y fx{t)=x 



(4.11) 



P{t,x;y\to,xo;yo) 



[S)u{t)][&f,m&xm&y{t)] 



'y{to)=yo Jx{to)=xo 
X exp 

X r (x(t) - fix{t)) - v{t))5^{m - h{x{t)) - ^(t)). 
The delta functional integrations over v{t) and /i(t) are trivial, so that 



(4.12) 



V ^"-^ i=l ^'^M i=i Jto ^ Jt^ OX^ ^ 



P{t,x;y\to,xo;yo) 



y{t)=y fx(t)=x 



y(to)=yo Jx{to)=xo 



i^xm^ymeM-s), 



where 



(4.13) 



(4.14) 



S is referred to as the action. 
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4.2 Sampled Continuous Measurements 

The path integral formula, Equation 4.13| , can be further simplified by considering the case 
when the measurements are available at discrete time instants. Specifically, suppose measure- 
ments are available at time and ti, and that the transition probability density at time 
ti is desired. Further, assume that there are no measurements available between and tj. 
Then, 



-W / - hk{x{t))f = -— j dtY, [yl{t) + hl{x{t)) - 2hk{x{t))yi{t)] . 

fJ' J ti—\ k = l ^ ^i — l k = l 

(4.15) 

The quantity of interest is the state and so all state-independent contributions may be ignored. 

As the first term is independent of the state variables, it is can be absorbed in the 
measure. The second term can be added to the action term that is independent of y{t). It 
remains to evaluate the third term: 

— / ' Y.^k{x{t))yk{t)dt. (4.16) 

There are two issues in this evaluation. Firstly, this can be evaluated via the usual integration 
by parts, but it is important to note that it is valid only for symmetric discretization. Secondly, 
since the measurements are sampled, there are two possible interpretations: 

^ /'■ t ".(^mnm - 1 1 ^'r- '''<^-(*'«'(')*- (4.17) 



fc=l 



titZ~XT=ihk{xmk{t)dt. 



This leads to two possibilities: 



1 \K:^k=ihk{^[xiti) + x{ti^i)])[yk{ti)-yk{ti-i)], 

EfcLi hk {2[x{ti) + x{ti^i)]) [yk{ti^i) - ykiU^i)]- 



(4.18) 

Therefore, the path integral formula simplifies to 

P{ti,Xi\yi\ti-\,Xi-\\yi-\) = P{ti,Xi\ti-i,Xi^i) (4.19) 

exp (-^ Efc=i ^fc {^[x{ti) + x{ti^i)]) [ykiU) - ykiU-i)]) , 
^ exp ( ^ Ylk=i {^[x{ti) + x{ti^i)]) [yk{ti-i) - yk{ti-2)] 

where 

P(ti,Xi\U_i,x,_i)= ' ' [&x{t)]eM-S{U.i,U)), (4.20) 

Jx{ti-l)=Xi-l 
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and 



1 

S{ti_i,ti) = - dt 



1=1 



1=1 



k=l 



. (4.21) 



It shall be shown in Section ^ that these lead to the two forms of the Yau algorithm. 

5. Path Integral Formula: Explicit Time Dependence 

In this section, the fohowing model is considered: 

j dx{t) = f{x{t),t)dt + e{t)du{t), x{0)=xo, Q(t) G M"""", 
\dy{t) =h{x{t),t)dt + n{t)dix{t), 2/(0) = 0, ii(t) e M™^"". 



(5.1) 



The discussion is similar to that in Section ^ However, imposing a delta function constraint 
requires a different method since e{t) and n(t) are not square matrices, and hence are not 
invertible. 

5.1 General Result 

As in Section ^, the transition probability density is computed by averaging over the signal 
and measurement ensembles, i.e., 

Pit,x;y\to,xo;yo) = j [dp{u{t))][dp{ix{t))] (5.2) 

X [dx{t)]J5''[x{t) - f{x{t),t) - e{t)u{t)] 
X [dmW'^im - h{x{t),t) - n{t)t,{t)] 
X 5"[x(t) - x]U(i„)=^.3 5™[y(t) - y]\y(t,)=y,- 

From the assumptions of the signal and measurement noise processes, it is evident that 

p 



[dp{v){t)] = [^Ki)]exp {Q-\t))^^v,{t)dt^ , 

[dp{p){t)] = [^M(t)]exp l^-i {W~\t))^pS)dt^ ■ 

The Jacobian J follows from the functional derivative of the Langevin equation: 



(5.3) 
(5.4) 



5xj{t') 



i{t) - fi{x{t),t) -Yeiait)Uait) 



a=l 



J. d dfi , , s . 



dt dx^ 



d_ 
dt' 



6{t-t'), 



(5.5) 



S,.8{t-t')-9{t-t')^{x{t),t) 
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Hence, 



J = det det [5{t - t') - e{t - t')^{x{t),t)^ , 

= J^det (^6{t - t') - e{t - t')^{x{t),t))^ , 



(5.6) 



where ryV is an irrelevant constant, or. 



In J = In det 



(5.7) 



The Jacobian Jy is trivial (as y(t) is y— independent) and can be absorbed into the measure. 
Thus, so far. 



P(t,x;y|to,a;o;2/o) 



yito)=yo Jx{to)=xo 



X 5" U 



n p 



[i^x(t)][^z.(t)]exp / E E ^amQab{t)r'n{t)dt 



* ii=l a,b=l 



{t) - fiix{t),t) -Y,^eia{t>a{t)^ exp l^-i I J2 
S"" (yk{i) - hk{x{t),t) - Y,^Ut)^^c{t) \ . 



dx> 



(5.8) 

\x{t),t)dt 

(5.9) 



c=l 



Using the Fourier integral version of the delta function (Equation |3^ 



Pit, x; y\to,xo; yo) = / [^x{t)] [&y{t)] [^u{t)] [^^(t)] (5.10) 



[^A(t)]exp I f^X^it) 



i{t) - fiix{t),t) -Y,eia{tHit) 



a=l 

q 



[^K(t)]exp i / ^Kk{t)[yk{t) - hk{x{t),t) - Y,'^kc{t)fic{t)] 



k=l 



c=l 



X exp 



II Yl Mt){Q'Ht))abn{t)dt] xexp 1-^ I E 

a,b=l I \ i=l 



dxi 



{x{t),t)dt 



Integrating over i^{t), and leads to 
ry('t)=y fx{t)=x 

P{t,x;y\to,xo;yo) 



[^A(t)] [&X][&k] exp 



v{to)=yo Jx{to)=XQ 



dt 



i=l 



(5.11) 



X exp 



dt 



dt 



5^ Ht)eia{t)Qab{t)el^{t)X,{t) - i^Xiit) {x^{t) - fi{x{t),t)) 

i,j=l a, 6=1 



j=l 



X exp 

_k,l=l c,d=l 

Integrating over Xi{t) and Kfc(t), it is clear that 

ry{t)=y f-x{t)=x 
P{t,x;y\to,xo;yo) = / 

■Jy(to)=yo Jx{to)=xo 

where the action is given by 

n p 



X] X] '^di)nkc{t)Qcd{t)n%{t)Ki{t) - i^Kkit) {yi{t) - hi{x{t),t)) 



k=l 



?x(i)][^y(t)]exp(-5), 



(5.12) 



-1 /" 

^=2^^^ ^ " ^^^^^^^^'^^^ (eia(t)Qa6(0eSW)"' " (5-13) 



+ 



1 



dt 



dfi 



(y^W - hk{x{t),t)) {nUt)Wcd{t)nl{t)) {yi{t) - /i;(x(t), f)) + ^ ^(x(t), t) 



i=l 



_k,l=lc,d=l 

5.2 Sampled Continuous Measurements 

As in Section an attempt is made to reduce the path integral to just the state variables. 
The general path integral formula (Equation 5.12| ) cannot be simplified in the sampled mea- 
surement case as easily as in the time dependent case unless some additional assumptions are 
made. 

First, observe that if e{t)Q{t)e^ {t) and n{t)W{t)n{t) are hylnxn and h^Imy^m and h{x{t)) 
is not explicitly time dependent, the resulting path integral is the same as Equation 4.19| . This 
is the case studied in 

Secondly, if the conditions above are relaxed to allowing explicit time dependence of the 
drift term in the state model, i.e., f{x{t),t), and and C = {n{t)W {t)n'^ (t)) , then the path 
integral formula becomes 

P{ti,Xi;yi\ti_i,Xi-i;yi^i) = P{ti,Xi\ti-i,Xi^i) (5.14) 

exp((^ EZ=i ^fc ihH^^) + ic-')ki - yk{U-2)]) , 



where 



?x(t)]exp(-5(t,_i,i,)), 



(5.15) 
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and the action is given by 



dt 



ti-1 



1 = 1 



1 = 1 



t "I 



^Y,{±,{t)-f,{x{t),t)f + Y,^{x{t),t) + ^ / dtY,hi{x{t)) 



■^l J to 



k=l 



(5.16) 



Finally, consider the case when e{t)Q{t)e{t) is time-independent and invertible, but oth- 
erwise arbitrary, and C = {n(t)W {t)n^ (t)) . Note that C is a constant, symmetric matrix and 



assumed to be invertible. The path integral formula is given by Equations 5.14 and 5.15 and 
with the action S'(tj_i,tj) is given by 



rti " P 

S{ti.i,ti) = 1^1 dtY^Y. (^'(*) - /*(^(*)'*)) {(^ia{t)Qab{t)elj{t)y^ {ij{t) - fjix{t),t)) 
■^ti-1 ^h^l 



(5.17) 



h,{x{t)) (C7-i),^ hi{x{t)) + ^.i^(t),t) 



k,l=l 



i=l 



In Section 6.2, the partial differential equation satisfied by P{t, x\to, xq) is derived. This is 
a straightforward generalization of the YYe to the state model with explicit time dependence. 
However, it is not rigorously shown that it approximates the solution of the DMZ equation 
for this case. 



6. Partial Differential Equations Satisfied by the Path Integral Formulas 
6.1 Time-Independent Case and the YYe 

In this section, it is proved that the path integral formula derived in the time-dependent case 
satisfies the YYe. Our discussion follows the method used by Feynman to verify the path 
integral formula he derived for the fundamental solution of the Schrodinger equation fP^ . 

From the Chapman-Kolmogorov semi-group property (or equivalently properties of the 
Gaussian measure), it follows that 

P{t + e, x\to, xo) = j P{t + e> ^l*> x')P{t, x'\tQ,XQ){d'^x'], (6.1) 

where P{t,x\tQ,XQ) is given by Equation 4.20| . 



Now according to Equation 4.20 



P{t + e,x\t,x')= (6.2) 

\ i=l i=l ^ k=l / 
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where x = ^{x + x'). The dominant contribution is when the following condition is satisfied: 

x-x' - ef{x) ^ 0. (6.3) 

In this region this may be written as 

X = x' + e/(x) + r], or (6-4) 
X = x' + ef{x) + r], 

where the equalities are valid to 0(e). Substituting this into Equation S^, leads to (keeping 
terms up to 0(e)) 

P{t + e,x\to,xo) = A f exp(--^^?7n (g^) 



2h^e . ^ 

4 = 1 



m 



i=i ""/^fe=i 



^ / exp 



\ i=l k=l / \ i=l / 

P{t,x - ef{x,t) - r]\to,xo). 

Note that the Jacobian of the the transformation from x' to rj to 0(e) is included in the 
second step. 

The constant A is fixed by 

A l^exp J-^r?2^ {d-r?} = 1. (6.6) 

Hence, it follows that 

A J r]irjj exp (^"^ E ^'^''^^ = ^'^^'^^i" (6-^) 
The left hand side of the equation is 

dP 

P{t, x'\to,xo) + e—{t, x\to, xo). (6.8) 



The second order expansion (in rj) of P{t, x'\tQ, xq) 
P{t, x'\to,xo) = P{t, X - ef{x, t) - r]\to, xq), 



(6.9) 



i=l * ^ k=l / 



dP 1 d P 

^(e/i(x) +r]i) —{t,x\to,xo) + 2 X] g^.g^ . ^I^O' ^o)- 



Only in terms of 0(e) are of interest here. The only terms nonvanishing are terms linear in 
e and quadratic in ij. Then it is easy to verify that the diffusion part of the Yau Eqution 
follows from the term quadratic in 77, so that 



dP 
'dt 



{t, x\to,xo) = ^ ^ 'fTT^^' x\to,xo) - ^ 



i=l 



d_ r 

dx 



fi{x)P{t,x\to,xo) 



1 

2h, 



J2hkix)Pit,x\to, 



k=l 



(6.10) 



K A d^P . 
2 ^ dxf 

P{to,x\to,xo) = 5"-{x - xo). 

Hence, P satisfies the YYe with a delta- function initial condition^. 
6.2 Time-Dependent Case 

We now derive the PDE satisfied by the path integral formula derived for the time-dependent 
case. This is the analog of the YYe in the time-dependent case. 

From the Chapman-Kolmogorov semi-group property, it follows that 



Xo) 



P{t + e,x\tQ,xo) = j P{t + e,x\t,x')P{t,x'\to,xo){d''x'}. 



.11 



According to Equation 5.14 and 5.17| 

P{t + e,x\t,x') = 
/ 1 " 

^exp I - ^ X] [x-x' - e/i(x,t)] [eia{t)Qab{t)elj] [x - x' - ef{x,t)] 
\ ^ i,j=i 

^ i=i ^> k,i=i J 

The contribution is negligible unless 

X — x — ef{x, t) « 0. 



.12) 



(6.13) 



^It is straightforward to see that P{to,x\to,xo) = lim^^o P{t + e,x\t,xo) = S"{x — xo). 
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so that 



X = x' + ef{x, t) + 77, or 
X = x' + ef{x, t) + 7]. 



(6.14) 



where the equahties are vahd to 0(e). Substituting this into Equation 3.2 , 

n p 

2e 



/•OO / 1 " ^ 

P(t + e,x|to,xo) = exp I -— ^ ^ m{eia{t)Qab{t)ebj{t)y^rij 



(6.15) 



i,j=l a,b=l 



1=1 



e 

2^: 



/ifc(x) (C-i),^/iKx) I P{t,x'\to,xo){d^x} 

k,l=l ^ 

A exp -- ^ ^ ??i(e^a(t)Qafe(i)ef,(t))-'??j 

\ i,j=la,b=l j 



P{t, x'\to, xo) ( 1 - I E |^(x, t) ) {d"ry} 



i=l 



j-oo I -n p \ 

^ / E E ^i(e*a(i)Qafe(i)ef,(t))-^7?, 

V M=la,fe=l / 



i=l 



hk{x) {C-')^^hi{x) P(t,x-e/(x,t)-7?) 



fci=i 



where the constant A is fixed by 



POO / 1 " ^ \ 

A / {d>}exp --EE ^.(e.a(t)Qa6(t)e^,(t))-ir?,- = 1, 



(6.16) 



so that 



{d"??} r/fcr/i exp I -^E E Viieiait)Qab{t)eljit))-^Vj 1 = ^ E {'^kait)Qab{t)eli 



i,j a, 6=1 



a,fe=l 



(6.17) 



The Jacobian of the the transformation from x' to rj to 0(e) is included in the second step in 
Equation S.15. 

The left hand side of the equation is 



dP 

P{t, x'\to, xo) + e-^(i, x\to, Xo). 



(6.18) 



To second order in rj, P{t, x'\tQ, xq) is 

P{t, x'\to, xo) = P{t, X - ef{x, t) - r]\to,xo), 

= Pit, x|to, xo) f 1 - e t) - ^ 



(6.19) 



k,l 



dP 1 " d^P 

^{efi{x,t)+r]i) —{t,x\to,xo) + ^i'nj q^,q^ i't^x\tQ,XQ). 



Since only terms of 0(e) are of interest, the relevant terms are linear in e and quadratic in 
r/.Thus, the PDE for the transition probability density is 



' dP ^ n p Q2p 



dt 



■{t,x\to,xo) 



i,j=l a,b=l 

d 



^ r _ 



i=l 



(6.20) 



k,i=i 

P{to,x\to,xo) = d^'ix - Xo). 



^ m 

-ijtY. ^fc(^) (C-')fc^/i/(x(t))P(t,x|to,xo), 



Clearly, this reduces to the YYe when e{t)Q{t)e^ (t) = I = hy = 1 and the drift in the 
signal model is not explicitly time dependent. 



7. Derivation of the Yau Algorithm 



In Section |3.3| it was noted that if i;j_i(tj_i, Xj_i) is the conditional probability density at 
time then the conditional probability density at time ti is given by 



Vi{ti,Xi) = / P{ti,Xi;yi\ti-i,Xi-i]yi-i)vi-i{ti_i,Xi-i) {(f xi-i} . 



(7.1) 



First consider the case discussed in Section H. Then 



Vi{ti,Xi) = (7.2) 

/ P{ti,Xi\ti-i,Xi-i) exp (-^ YJk=i {\[x{ti) + 2;(tj_i)]) [ykik) - yfc(ti-i)]) Wi-i(ti-i, Xj_i) {(Pxi-i] 
J P{ti,Xi\ti-i,Xi-i) exp ^--i Y^Z=i {^[x{ti) + x{ti^i)]) [yk{ti-i) - yk{ti^2)\^ i;(tj_i, Xi_i) {d"-Xi^i} . 



When ti — ti-i is small 



hk Q[x(ti) + x(t,_i)]^ ^ hk{x{ti^i)) + 0{^e), 



(7.3) 



- 30 - 



and Equation TJl becomes 



/ P{ti, Xi\ti^i,Xi-i) exp (-^^YJk=i hk{x{ti-i))[yk{ti) - yk{ti-i)yj Xj.i) {(Txi^i) 

J P{ti,Xi\ti^i,Xi-i) exp (^-^ Y^k=i hk{x{ti-i))[yk{ti-i) - yk{ti-2)]^ Vi^i{ti-i,Xi^i) {(Pxi^i} 

(7.4) 



In Section ^, it was already proved that is the fundamental solution of the 

YYe. This implies that Vi{ti,x) is the solution at ti of 



dt 



2 ^ dxj 

1=1 ' 



dxi 



1=1 * \i=l 

exp (yyjLi ^ji^) [yjiU) - VjiU-i)]^ ui{ti,x), 
exp (EjLi hj{x) [yjiti^i) - yj{ti^2)]) ni(Ti,x), 



1 ™ \ 
+ ^Ylhf{x)j Vi{t,x), 



Vi{ti-l,x) = Vi-l{ti-l,x) 



(7.5) 



This is precisely the Yau algorithm. 

Likewise, it is straightforward to see that for the general case studied in Section |5|, the 
Yau algorithm is extended to this case as follows: Vi{ti, Xi) is the solution at ti of the PDE 



dvi 
dt 



n p 



i,j=l a, 6=1 

d 



dxidxj 



{t,x) 



E Q^. [Mx^t)Vi{t,x)] 



i=l 



^ m 

— hk{x) {C~')j^^hi{x{t))v,{t,: 



(7.6) 



k,l=l 



Vi{ti^l,x) 



exp {J2]^k=if^ji^) (C ^)jk iVkiti) - yk{ti~i)]j Vi-i{ti^i,x), 
exp ( E™fc=i^i(a;) (C""^)jfc [yk{ti-i) - ykiti-2)]) Wj-i(tj-i,x). 



8. Some Remarks 



Following are come remarks on the Feynman path integral solution of the filtering problem: 

• Note that the Feynman path integral formulation has given a complete and self-contained 
solution of the continuous-continuous filtering problem. For instance, the DMZ equa- 
tion or its variants were not used as an input. On the contrary, the Feynman path 
integral formula naturally led us to the YYe and the Yau algorithm. 

• Unlike the YYe, the PI formula is valid even for time-dependent case. Thus, in principle, 
one can compute the conditional transition probablity density using the conventional 
methods (see, for instance, |p3||). 



- 31 - 



The YYe can be viewed as a local expression of the path integral formula. That is, a 
path integral is a global object, while the PDE is a local one. 

In this paper, the signal noise and measurement noise are assumed to be additive. This 



is a stronger condition than the orthogonalilty of the diffusion vielbein assumed in |13]. 
In a subsequent paper, this case will also be covered. 

It is clear that the robust DMZ equation can also be solved using the transition prob- 
ability density. However, the quantity of interest, the conditional probability density 
a(t,x) is computed directly by the path integral formula. 

It is noted that other algorithms can also be solved using the Feynman path integral 
formulas with obvious changes. For example, the splitting-up approach requires the 
solution of the FPKfe, which corresponds to the case h{x) = 0. Note that the FPKfe 
arises naturally in the solution of the continuous-discrete filtering problem |4C]. Of 
course, as noted previously, it would be unnecessary since the Yau algorithm has the 
best numerical properties. What is interesting to note is that the path integral formula 
naturally leads to the algorithm with the best numerical properties. 



9. Examples 

9.1 The Dirac- Feynman Approximation 

From the discussion in the previous sections, it is evident that computing the path integral 
requires computing the Lagrangian, L, defined by 5 = J dtL{x,x,t). When the time step is 
infinitesimal, it is given by Equation 90. Therefore, the simplest and crudest approximation 
is to use this for non-infinitesimal time step. Specifically, the fundamental solution of the 
YYe is approximated as follows: 



P{t",x"\t',x') : 

{t" -t') \ 1 



exp 



1=1 



I I r 

II _ -j-l J'^ 



x" + x' 



2 n 



^ dxi 



i=l 



X" + x' 



^ m 



k=l 



(9.1) 

x" + x'" 



9.2 Example 1 

As an example, consider the following continuous-continuous filtering model that has been 



(9.2) 



studied in Q] (and |8|) 

dx{t) = acos{bx{t))dt + axd\i{t), 
dy{t) = tan~^ x{t) + ayd\M{t). 

The Lagrangian for this model is given by 

— 7T (x(t) — acos(bx(t)))'^ sm(bx(t)) -\ ^ (tan^"*^ x(t))' 

2az 2 2af, 



(9.3) 
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The parameters chosen were as in the reference. Specifically, a = 1.2, b 
grid spacing Ax = 0.01 and extent [—1.5, 1.5], temporal grid spacing At 
steps. 



= 3,crx = 0.3, spatial 
= 0.01 with 200 time 



A sample of the signal and measurement state processes with o,(=0.3, and 0^=0.05. 




100 120 



180 200 



Figure 1: A sample of state and measurement processes given by Equation 9.2 



In the first set, the measurement noise was set as ay = 0.05. Figm'e |l| shows a sample of 
state and measurement processes. In Figure ^ is plotted the conditional mean computed using 
Equation |9.1| in the Yau algorithm^. Also plotted are 2a limits. The conditional mean and 
variance were computed from the computed conditional probability density. The fact that 
the target was mostly within the 2a limits of the conditional mean shows that the tracking 
performance of this algorithm is quite good. 

In the next set, the measurement noise was set as o"„ = 0.0125. For this "small noise" 



case, most of the algorithms studied in |44| failed. In Figure |^ is plotted a sample of signal 
and measurement processes. The conditional mean and the 2a limits computed using the 
Yau algorithm for this instance is plotted in Figure 

It is seen that good tracking performance is maintained for the small noise case even 
when the crudest path integral approximation is used in the Yau algorithm. 



Since there was negligible difference in performance between the pre-measurement and post-measurement 
forms, only the former was employed. 
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Conditional Mean with 2o bounds for the sample path (measurement noise 0^=0.05). 



1.5 



0.5 



-0.5 



-1 



-1.5 




3 ^ 



20 40 60 80 100 120 140 160 180 200 



Figure 2: Conditional mean and 2a limits computed using the Yau algorithm with path integral 
approximation of Equation 



9.3 Example 2: Cubic Sensor Problem 



The cubic sensor problem is defined by tlie following signal and measurement model: 



d\{t) = axd\i{t), 

dy{t) = x^{t)dt + cjyd\M{t). 



(9.4) 



It is a well-studied nonlinear filtering problem because it is one of the simplest example 
examples of a filtering problem that is not finite dimensional (see, for instance, [^] and 
references therein). 



For the simulation of the cubic sensor problem, the following model parameters were 



- 34 - 



A sample of the signal and measurement state processes with Oj,=0.3, and 0^=0.0125. 




200 



Figure 3: A sample of state and measurement processes given by Equation |9.2| — small measurement 



noise case. 



chosen (as in |^ 



At = 0.001, 
Ax = 0.01, 
xo = iV(0,0.01), 

CTr = 1, 



ay = 1, 



(9.5) 



T = [0,20]. 

The Extended Kalman filter, a sub-optimal filter which approximates the conditional 
probability by a Gaussian, for the cubic sensor problem is given by 

dx = 3x^P{dy - x^dt), 
dP = {l- 9x'^P^)dt. 

The Lagrangian for the cubic sensor problem is 



(9.6) 



(9.7) 
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Conditional Mean with 2o bounds for the sample path (measurement noise o =0.01 25). 

1.5 



1 



0.5 







-0.5 



-1 



-1.5 

20 40 60 80 100 120 140 160 180 200 



Figure 4: Conditional mean and 2cr limits computed using the Yau algorithm with path integral 
approximation of Equation with small measurement noise. 

The simplest approximation is to use the simplest approximation of the Lagrangian 

= exp (-^ [^)'- ^ . ,9.8, 

Figure ^ shows the performance of the approximate path integral. Specifically, the con- 
ditional mean along with the one standard deviation bounds computed using the computed 
conditional probability density is plotted. Observe that the EKF fails completely in this case. 
As noted in |^5| , this is because the EKF considers only the first two moments (which vanish 
here) ; it is the fourth central moment that plays a crucial role in this example (for the chosen 
initial condition). Also, note that the state is within the la region for most of the time. This 
shows that, unlike the EKF, the path integral filter has a reliable error analysis. 

After an initial period, the performance of the path integral approximation is seen to be 
excellent and comparable to that obtained using PDF methods in |^5[. However, the crucial 
point is that the path integral solution is equally simple for the higher- dimensional case with 
more complicated models (e.g., colored noise), whereas a PDE solution would he significantly 
harder, if not impossible, to implement in real-time. 
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Conditional Mean with 1 a bounds for thie sample path (measurement noise o =1 ). 




Figure 5: A simulation result of cubic sensor problem. 



9.4 Comments 



It is remarkable to note the very good performance is obtained using the crudest approxi- 
mation. Of course, when the time step is large, it will fail (unlike the path integral formula 
itself). However, the practical situation is that the time steps are often small. Therefore, this 
may not be a bad approximation in many cases. 

The implementation of this method is trivial. The contrast with other methods, such 
as those studied in |44], is striking. For instance, many of those methods require off-line 



computation of complicated partial differential equations with uncertain numerical properties. 

The results obtained in this paper used single time-step. More accuracy can be obtained 
quite simply using multiple time steps. Also, the computation of the transition probability 
density can be done off-line, but the on-line computation was not an onerous burden for this 
case. 

It is important to note also that the transition probability density matrix (or tensor in 
the general case) is sparse (determined by h^,hi,). This is of great importance in higher 
dimensional filtering problems because 

• Sparse matrix storage requirements are small, 

• The relevant transition probability density matrix elements can be computed based on 
the conditional density in the previous step, and 

• Sparse matrix computations are very fast. 
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Note that unlike some other approximation techniques studied in |^^, the conditional 
probability density is obviously always positive (provided, of course, that the initial distribu- 
tion is positive). 



10. Conclusion 

In this paper, the formal path integral solution to the continuous-continuous nonlinear filtering 
problem has been presented. The solution is universal, i.e., the initial distribution may be 
arbitrary. We verified that the path integral formula satisfies the YYe when the drift is not 
explicitly time-independent. Since the path integral measure is manifestly positive, positivity 
is maintained if the initial distribution is positive. 

A path integral formulation has several advantages. Path integrals have led to theoretical 
insights in other areas including quantum mechanics, quantum field theory and even math- 
ematics. It is demonstrated in this paper that it is quite easy to express the fundamental 
solution of the YYe in terms of path integrals. In fact, the YYe has been extended to the 
case when the drift is explicitly time-dependent. In a future paper, this analysis is extended 



to the case of multiplicative noise [ 46 1 



Finally, path integrals are very suitable for numerical implementation. Practical path 
integral filtering techniques will be presented in subsequent papers. 
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